import numpy
import scipy.sparse as sparse
import scipy.sparse.linalg as linalg
import matplotlib.pyplot as plt
import plotly.graph_objects as go
import pandas as pd
from ipywidgets import interact, SelectionSlider
import plotly.io as pio
import plotly.offline as py
pio.renderers.default = "notebook_connected"
py.offline.init_notebook_mode()

Numerical Analysis of 2D Diffusion-Reaction Equations Using an Implicit-Explicit Method: Solution and Visualization

• Overview: This is a showcase of a hybrid implicit-explicit method to solve the diffusion-reaction PDE in 2-dimensions. It uses Plotly to provide a stroboscopic and interactive visualization of the solution.

• Connection to data science: The PDE has indirect connections to time-series analysis and forecasting. The diffusion equation can be used to smooth time-series data, removing short-term fluctuations while preserving long-term trends. This is analogous to applying a low-pass filter. The diffusion helps in removing short-term fluctuations, and the reaction term can introduce additional effects such as trends or external influences.

• Why use a hybrid method?: The diffusion term involves second spatial derivatives, which tend to amplify numerical errors in explicit methods like Forward Euler. The stability criterion for the diffusion term is $\Delta t=\frac{\Delta x^2}{4D}$, where $D$ is the diffusion coefficient. This requires small time steps, making the method not time-efficient. By using an implicit method for the diffusion term and an explicit method for the reaction term, the hybrid approach is more time-efficient. ​

The Implicit-Explicit Method:

The PDE is

\begin{aligned} u_t &= \sigma D_1 (\partial_{xx}+\partial_{yy})u + f(u, v) \\ v_t &= \sigma D_2 (\partial_{xx}+\partial_{yy})v + g(u, v) \end{aligned}

with the reaction terms

\begin{aligned} f(u,v) &= \alpha u (1 - \tau_1 v^2) + v (1 - \tau_2 u) \\ g(u,v) &= \beta v + \alpha \tau_1 u v^2 + u (\gamma + \tau_2 v). \end{aligned}

The following numerical method solves the PDE over the square domain $\{(x,y):-1\leq x,y\leq1\}$ with periodic boundary conditions.

Numerical method:

\begin{aligned} U^\ast &= U^n + \Delta t \sigma D_1 (\partial_{xx}+\partial_{yy}) U^\ast \\ V^\ast &= V^n + \Delta t \sigma D_2 (\partial_{xx}+\partial_{yy}) V^\ast \\ U^{n+1} &= U^\ast + \Delta t f(U^\ast, V^\ast) \\ V^{n+1} &= V^\ast + \Delta t g(U^\ast, V^\ast). \end{aligned}

The following function distretizes the spacial derivative using a five-point stencil in 2d and returns a sparse matrix reprsentation.

def laplacian_discretization(m):
    """Constructs a sparse matrix that discretizes the 2d Laplacian
    Uses a five-point stencil and periodic boundary conditions.
    """
    delta_x = 2.0 / (m + 1)
    
    # Primary discretization
    e = numpy.ones(m)
    T = sparse.spdiags([e, -4.0 * e, e], [-1, 0, 1], m, m)
    S = sparse.spdiags([e, e], [-1, 1], m, m)
    I = sparse.eye(m)
    A = sparse.kron(I, T) + sparse.kron(S, I)
    
    # Construct periodic BCs
    e = numpy.ones(m**2)
    A_periodic = sparse.spdiags([e, e],[m - m**2, m**2 - m], m**2, m**2).tolil()
    
    # Left & right BCs:
    for i in range(m):
        A_periodic[i * m, (i + 1) * m - 1] = 1.0
        A_periodic[(i + 1) * m - 1, i * m] = 1.0
    
    # Combine two matrices
    A = A + A_periodic
    A /= delta_x**2
    A = A.todia()
    
    return A

def f_reaction(U, V, sigma, tau_1, tau_2, alpha, beta, gamma):
    return alpha * U * (1.0 - tau_1 * V**2) + V * (1.0 - tau_2 * U)

def g_reaction(U, V, sigma, tau_1, tau_2, alpha, beta, gamma):
    return beta * V * (1.0 + alpha * tau_1 / beta * U * V) + U * (gamma + tau_2 * V)

def forward_euler_step(U, V, delta_t, A, sigma, f, g, D1=0.5, D2=1.0):
    """Take a single forward Euler step on the reaction-diffusion equation"""
    
    U_new = U + delta_t * (sigma * D1 * A * U + f(U, V))
    V_new = V + delta_t * (sigma * D2 * A * V + g(U, V))
    
    return U_new, V_new


def backward_euler_diffusion_step(U, V, A, delta_t, sigma, D_1, D_2):
    """Take a single backward Euler step on the reaction-diffusion equation"""
    U = linalg.spsolve((sparse.eye(A.shape[0]) - delta_t * sigma * D_1 * A), U)
    V = linalg.spsolve((sparse.eye(A.shape[0]) - delta_t * sigma * D_2 * A), V)
    return U, V

def imex_solver(sigma, tau_1, tau_2, alpha, beta, gamma, D_1, D_2):
    # Alias reaction functions with the above parameters
    f = lambda U, V: f_reaction(U, V, sigma, tau_1, tau_2, alpha, beta, gamma)
    g = lambda U, V: g_reaction(U, V, sigma, tau_1, tau_2, alpha, beta, gamma)

    # Set up grid
    m = 150
    delta_x = 2.0 / m
    x = numpy.linspace(-1.0, 1.0, m)
    y = numpy.linspace(-1.0, 1.0, m)
    Y, X = numpy.meshgrid(y, x)

    # Initial data
    U = numpy.random.randn(m, m) / 2.0
    V = numpy.random.randn(m, m) / 2.0
    #fig = plt.figure()
    #axes = fig.add_subplot(1, 1, 1, aspect='equal')
    #plot = axes.pcolor(x, y, U, cmap=plt.get_cmap("viridis"))
    #fig.colorbar(plot)

    # Setup spatial discretization
    U = U.reshape(-1)
    V = V.reshape(-1)
    A = laplacian_discretization(m)

    # Time
    t = 0.0
    t_final = 30.0
    delta_t = delta_x / (10.0 * sigma)
    num_steps = int(numpy.round(t_final / delta_t))

    data = []
    data.append([x,y,U, t])
    
    # Evolve in time
    next_output_time = 0.0
    for j in range(num_steps):
        U, V = backward_euler_diffusion_step(U, V, A, delta_t, sigma, D_1, D_2)
        U, V = forward_euler_step(U, V, delta_t, A, sigma, f, g)   
        t += delta_t

        if t >= next_output_time:
            next_output_time += 5.0
            U_output = U.reshape((m, m))

            #fig = plt.figure()
            #axes = fig.add_subplot(1, 1, 1, aspect='equal')
            #plot = axes.pcolor(x, y, U_output, cmap=plt.get_cmap("viridis"))
            #fig.colorbar(plot)
            #axes.set_title("t = %s" % t)
            data.append([x,y,U_output, t])

    #plt.show()
    return data

# Parameters
data = imex_solver(sigma=0.0021, tau_1=3.5, tau_2=0, alpha=0.899, beta=-0.91, gamma=-0.899, D_1=0.5, D_2=1.0)
print('checkpoint')

checkpoint

# define time_intervals
time_intervals = []
for i in range(len(data)):
    time_intervals.append(data[i][3])


def extract_data(time):
    j = time_intervals.index(time)
    z_data = data[j][2]
 
    return z_data

def update_plot(time):
    
    z_data = extract_data(time)
    
    z_data = z_data.reshape((150, 150))
    

    fig = go.Figure(data=[go.Surface(z=z_data)])

    
    fig.update_traces(contours_z=dict(show=True, usecolormap=True,
                                      highlightcolor="limegreen", project_z=True),)
    fig.update_layout(title='Topographical 3-D Surface Plot of the Solution with Contours at t={}'.format(time), autosize=False,
                      scene_camera_eye=dict(x=1.7, y=1.3, z=0.64),
                      width=1000, height=700,
                      margin=dict(l=85, r=50, b=65, t=90),
    )

    fig.show()

for i in time_intervals:
    update_plot(i)
    
# These widgets are currently configured for .ipynb environments
# time_slider = SelectionSlider(options=time_intervals, value=time_intervals[5], description='Time')
# interact(update_plot, time=time_slider);

References

[1] Mandli, K. T. (n.d.). Reaction-Diffusion Demo. In Numerical Methods for Partial Differential Equations [Jupyter Notebook]. GitHub. Retrieved July 27, 2024, from https://github.com/mandli/numerical-methods-pdes/blob/master/reaction-diffusion_demo.ipynb

[2] Ketcheson, D. (n.d.). Finite Difference Course. GitHub. Retrieved July 27, 2024, from https://github.com/ketch/finite-difference-course

[3] LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations: Steady State and Time Dependent Problems. Society for Industrial and Applied Mathematics (SIAM).